Section: New Results

Computational Statistical Physics

Participants : Matthew Dobson, Claude Le Bris, Frédéric Legoll, Tony Lelièvre, Francis Nier, Stefano Olla, Grigorios Pavliotis, Giovanni Samaey, Gabriel Stoltz.

The extremely broad field of Molecular dynamics (MD) is a domain where the MICMAC project-team, originally more involved in the quantum chemistry side, has invested a lot of efforts in the recent years. Molecular dynamics may also be termed computational statistical physics since the main aim is to numerically estimate average properties of materials as given by the laws of statistical physics. The project-team studies both deterministic and probabilistic techniques used in the field.

Free Energy calculations

For large molecular systems, the information of the whole configuration space may be summarized in a few coordinates of interest, called reaction coordinates. An important problem in chemistry or biology is to compute the effective energy felt by those reaction coordinates, called free energy.

In the article [51] , Tony Lelièvre, Mathias Rousset and Gabriel Stoltz study the application of constrained Langevin dynamics to the computation of free energy differences, by thermodynamic integration techniques and fluctuation relation (à la Jarzynski).

The work by T. Lelièvre and K. Minoukadeh on the longtime convergence of the ABF method in a particular bi-channel scenario (which was already mentioned in last year's activity report) has been accepted for publication [50] . Likewise, the work by Nicolas Chopin (CREST, ENSAE), T. Lelièvre and G. Stoltz on application of the ABF method to Bayesian inference is about to appear, see [30] .

Sampling trajectories

There exist a lot of methods to sample efficiently Boltzmann-Gibbs distributions. The situation is much more intricate as far as the sampling of trajectories (and especially metastable trajectories) is concerned.

In [32] , T. Lelièvre and D. Pommier, in collaboration with F. Cérou and A. Guyader (INRIA Rennes, ASPI) investigated the interest of an Adaptive Multilevel Splitting algorithm to compute reactive paths, and estimate transition rates. The obtained results are very interesting. Current research aims at testing the technique on practical cases.

In [49] , C. Le Bris and T. Lelièvre together with M. Luskin and D. Perez have proposed a mathematical analysis of the parallel replica algorithm,introduced by A. Voter in 1997 to efficiently simulate metastable trajectories. This work opens a lot of perspectives using a generic tool (the quasi stationary distribution) to relate continuous state space dynamics (Langevin type dynamics) to discrete state space dynamics (kinetic Monte Carlo type models). A follow-up work consists in theoretically investigating another related approach, the hyperdynamics method.

Nonequilibrium systems

The efficient simulation of molecular systems is known to be a much more complicated problem when the system is subjected to a non-conservative external forcing than when the system experiences conservative forces. Together with the sampling of metastable dynamics mentioned above, these are the two major research focus in MD of the project-team.

On this topic, G. Stoltz continued his long lasting collaboration with physicists at CEA/DAM on reduced models for shock and detonation waves. More precisely, he published two works applying some simulation techniques he devised to actual energetic materials of interest to physicists, namely (i) a technique to sample constraints in average and allowing to compute the Hugoniot curve efficiently, which was applied to reacting TATB [24] ; and (ii) a reduced stochastic dynamics to model detonation waves, applied to a material with properties close to nitromethane and allowing an atomistic simulation of the shock-to-detonation transition [55] .

F. Legoll and G. Stoltz pursued their studies of the anomalous thermal conductivity of one dimensional chains. They have investigated the case of a chain of rotors subjected to a mechanical forcing. In collaboration with A. Iacobucci and S. Olla (CEREMADE, Paris Dauphine) they have shown in [44] that the mechanical forcing can have a counter-intuitive effect and reduce the thermal current. Besides this numerical study, G. Stoltz, in collaboration with C. Bernardin (ENS Lyon) considered in [18] the issue of thermal transport in one of the simplest possible one dimensional model, a chain of oscillators whose kinetic and potential energy functions are the same, and which are subjected to a stochastic noise exchanging all the variables. The system therefore has only two conservation laws, the energy and the total length. A hydrodynamic limit consisting of a system of conservation laws can be obtained before the onset of shocks. However, the thermal transport is anomalous: this can be proved by analytical computations for harmonic interactions, or demonstrated numerically in the general case.

G. Stoltz also studied techniques to compute the viscosity of fluids using steady state nonequilibrium dynamics with an external nongradient bulk forcing, in the framework of the context of the PhD of Rémi Joubaud, see [45] . The two authors have proved a linear response result, and obtained asymptotic scalings of the viscosity in terms of the friction coefficients of the underlying Langevin dynamics. G. Stoltz and G. Pavliotis are now extending the results to the case of time dependent nongradient external forcings.

Nonequilibrium molecular dynamics simulations can also be used to compute the constitutive relation between the strain rate and stress tensor in complex fluids. This is fulfilled simulating molecular systems subject to a steady, non-zero macroscopic flow at a given temperature. Starting from a bath model, M. Dobson, F. Legoll, T. Lelievre, and G. Stoltz have derived a Langevin-type dynamics for a heavy particle in a non-zero background flow.The resulting dynamics, which is theoretically obtained when a unique large particle is considered, is numerically observed to also perform well when a system of many interacting particles within shear flow is considered.

Effective dynamics

For a given molecular system, and a given reaction coordinate $\xi :{ℝ}^n\mapsto ℝ$, the free energy completely describes the statistics of $\xi \left(X\right)$ when $X\in {ℝ}^n$ is distributed according to the Gibbs measure. On the other hand, obtaining a correct description of the dynamics along $\xi$ is complicated.

F. Legoll and T. Lelièvre have continued their work on the definition and the analysis of a coarse-grained dynamics that approximates $\xi \left(X_t\right)$, when the state of the system $X_t$ evolves according to the overdamped Langevin equation (which is ergodic for the Gibbs measure). The aim is to get a coarse-grained description giving access to some dynamical quantities, such as residence times in metastable basins. These basins are usually assumed to be completely described by $\xi$. They have proposed an effective dynamics, which is derived using conditional expectations. The first accuracy result, obtained in 2010, is phrased in terms of an estimate on the relative entropy between the law of $\xi \left(X_t\right)$ and the law of its approximation, at any time $t$ (see [63] for a description of the results in a simple case). If an appropriate time-scale separation is present in the system, then the effective dynamics is accurate in the sense of time-marginals. The obtained numerical results showed that this dynamics can also be used to accurately compute residence times in potential energy wells, and thus seem to be accurate in a much stronger sense. Together with S. Olla, they have started to analyze the pathwise accuracy of the proposed coarse-grained dynamics.

The extension of the numerical strategy to the case when the reference dynamics on the whole system is the Langevin dynamics has also been studied. Promising numerical results have already been obtained in collaboration with a master's student (F. Galante).

Convergence to equilibrium

An important question for the analysis of sampling techniques is the rate of convergence to equilibrium for stochastic trajectories.

F. Nier continues to investigate the spectral properties of Witten Laplacians and Kramers-Fokker-Planck operators. In a recent collaboration with D. Le Peutrec and C. Viterbo, the Arrhenius law for metastable states, in its refined version also known as Eyring-Kramers law, was extended to $p$-forms. Very accurate analytical results have been provided for the exponentially small eigenvalues of Witten Laplacians acting on $p$-forms.

F. Nier has started a collaboration with T. Lelièvre about accurate exit laws for Smoluchowski processes, via a Witten complex approach. Investigations with G. Pavliotis and T. Lelièvre have started about non gradient diffusion systems.

Hamiltonian dynamics

Constant energy averages are often computed as long time limits of time averages along a typical trajectory of the Hamiltonian dynamics. One difficulty of such a computation is the presence of several time scales in the dynamics: the frequencies of some motions are very high (e.g. for the atomistic bond vibrations), while those of other motions are much smaller. This problem has been addressed in a two-fold manner.

Fast phenomena are often only relevant through their mean effect on the slow phenomena, and their precise description is not needed. Consequently, there is a need for time integration algorithms that take into account these fast phenomena only in an averaged way, and for which the time step is not restricted by the highest frequencies. In [38] , M. Dobson, C. Le Bris, and F. Legoll developed integrators for Hamiltonian systems with high frequencies. The integrators were derived using homogenization techniques applied to the Hamilton-Jacobi PDE associated to the Hamiltonian ODE. This work extends previous works of the team. The proposed algorithms can now handle the case when the (unique) fast frequency depends on the slow degrees of freedom, or when there are several fast constant frequencies.

Another track to simulate the system for longer times is to resort to parallel computations. An algorithm in that vein is the parareal in time algorithm. It is based on a decomposition of the time interval into subintervals, and on a predictor-corrector strategy, where the propagations over each subinterval for the corrector stage are concurrently performed on the processors. Using a symmetrization procedure and/or a (possibly also symmetric) projection step, C. Le Bris and F. Legoll, in collaboration with X. Dai and Y. Maday, have introduced several variants of the original plain parareal in time algorithm [35] . These variants, compatible with the geometric structure of the exact dynamics, are better adapted to the Hamiltonian context.